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Abstract 

The large time dynamics of a periodically driven Fokker-Planck process possess- 
ing several metastable states is investigated. At weak noise transitions between 
the metastable states are rare. Their dynamics then represent a discrete Marko- 
vian process characterized by time dependent rates. Apart from the occupation 
probabilities, so-called specific probability densities and localizing functions can 
be associated to each metastable state. Together, these three sets of functions 
uniquely characterize the large time dynamics of the conditional probability 
density of the original process. Exact equations of motion are formulated for 
these three sets of functions and strategies are discussed how to solve them. 
These methods are illustrated and their usefulness is demonstrated by means 
of the example of a bistable Brownian oscillator within a large range of driving 
frequencies from the slow semiadiabatic to the fast driving regime. 
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1. Introduction 

Chemical reactions have provided ubiquitous and versatile examples of ac- 
tivated transitions between two metastable states, formed by the reactants and 
products. In a chemical reaction the energy necessary for the activation most 
often stems from the (classical or even quantum mechanical) thermal energy 
that may accumulate in a single reaction coordinate and finally enable a tran- 
sition from reactants to products [1, 2, 3, 4]. In contrast to these thermally 
assisted escape processes other additional sources of energy may externally be 
provided for example by driving a system with metastable states by periodic 
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forces. Such periodically driven stochastic systems present a particular class 
of nonequilibrium processes that exhibit a broad variety of fascinating effects 
[5, 6, 7] such as stochastic resonance [8], directed transport of Brownian parti- 
cles in ratchet type periodic potentials [9, 10, 11] or other anomalous transport 
properties as for example negative mobility [12]. Apart from an external pe- 
riodic driving, these systems typically are subject to nonlinear dynamical laws 
and additionally experience fluctuating forces describing the random impact of 
the environment of the considered system [13]. Without the fluctuating forces 
the presence of nonlinearities often renders these systems multistable, i.e. such 
systems may approach different attractors [14], depending on their initial states. 
In combination with weak fluctuating forces these attractors become metastable 
states, which means that the system will be found most of the time in or close 
to one of these states while transitions between these states present rare events. 

Each of the principal constituents of the dynamics of a periodically driven 
nonlinear stochastic system is characterized by typical time scales such as the 
correlation time of the fast random forces (ff), Tfj, relaxation times r of the 
deterministic part of the dynamics, the period T of the driving force and the 
times Tms of typical sojourn within the different metastable states (ms). In 
this work we will assume that the correlation times of the fluctuating forces are 
much shorter than all other time scales such that a Markovian description of the 
dynamics is appropriate. Hence, we model the fluctuating forces by white noise 
(rff = 0) which moreover will be assumed to be Gaussian and weak. As a conse- 
quence of these assumptions the characteristic sojourn times of the metastable 
states are finite but much larger than any of the deterministic characteristic 
times (Tnis ^> t) [3]. This time scale separation implies that the transitions 
between the metastable states constitute a discrete Markovian process which 
will be investigated in more detail in the present work. We will demonstrate 
that this discrete process forms the backbone of the original continuous process 
on time scales that are much larger than the deterministic relaxation times r . 

Finally, the magnitude of the driving period T in relation to the determin- 
istic time scales r has a decisive influence on the system's dynamics. In the 
so-called semiadiabatic limit [15] the driving period is large compared to typical 
deterministic relaxation times independently of how large the driving period is 
compared with the typical sojourn times. Then the time-dependent transition 
rates are given by the frozen rates, i.e. their time dependence only results from 
the slow change of those system parameters that are varied by the driving pro- 
cess [16]. Within this framework stochastic resonance [17] and the dynamics of 
neuron models [18] have successfully been described. 

Outside the regime of the so called semiadiabatic limit the escape rates no 
longer instantly follow but rather lack behind the periodic driving [19]. In the 
present paper we investigate this regime of intermediate to fast driving in more 
detail and present effective methods to characterize the large time behavior of 
periodically driven Fokker-Planck processes with metastable states. 

Previous works on periodically driven processes with metastable states most 
often have been focussed on particular aspects such as on the dependence of the 
average life time of a metastable state [20, 21], of the exponentially leading part 
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of escape rates within linear response theory [22], or on rates in the weak noise 
Umit [23, 24]. 

We close this Introduction with a short outline of the paper. In Section 2 

wc introduce some important concepts of the deterministic dynamics of a peri- 
odically driven system with coexisting attractors. In Section 3 two alternative 
formulations of the conditional probability density function are presented for 
events that arc separated by a time that is much larger than the characteris- 
tic deterministic time r. The first form originates from the Floquet represen- 
tation of the conditional probability density of a periodically driven Markov 
process [5, 6] while the second expression explicitly refers to the dynamics of 
the metastable states. This second expression in particular contains quantities 
that characterize specific probability densities for each metastable state as well 
as localizing fimctions that allocate probabilities to the metastable states given 
the state of the full continuous system. In Section 4 we find equations of motion 
both for these metastable state specific probability densities and the localizing 
functions by comparing the two formulations of the conditional probability den- 
sity at large times. In Section 5 the theory is exemplified and numerically tested 
for a bistable Brownian oscillator. Section 6 closes with a summary. 

2. Characterization of the deterministic dynamics 

In the deterministic limit the considered system is described by the motion 
of a state x in a d dimensional state space S governed by a set of d coupled 
differential equations 

x = f(x,i), (1) 

where the vector field f(x, t) periodically depends on time with period T, i.e. 
f{x,t + T) = f(x, i). We denote the trajectory emanating at the time s from 
the point y by X(<|y, s) and assume that in the asymptotic limit of large times 
the motion is bounded and characterized by a set of n > 2 different attractors 
Aa{t) C S, a = 1 . . . n, such that each trajectory approaches either of the attrac- 
tors depending on its initial state and starting time, i.e. X(i|y, s) ^ x e Aa{t) 
for f — s sufficiently large. This relaxation process happens on a characteristic 
deterministic time scale of the considered system. The attractors periodically 
depend on time, i.e. 

Aa{t + T)=Aa{t) . (2) 

To each attractor a domain of attraction I'a(s) exists that consists of all states 
y at time s from which the a"^ attractor is reached. It is formally defined as 
T^a{s) = {y|X(t|y, s) £ Aa{t) for t — s — > oo}. At each fixed time the domains 
of attraction form a partition of the state space into disjoint subsets, which in 
general periodically depend on time 

Vc,{t + T)=Va{t). (3) 
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3. Conditional probability density of time-periodic Fokker-Planck pro- 
cesses with metastable states 

S.l. Floquet representation 

In many cases the description of a system in terms of deterministic equations 
of motion is sufficient in order to determine the typical behavior of the system 
with sufficient accuracy. However, the presence of weak random perturbations, 
which often can be modeled by Gaussian white noise, causes different effects 
depending on the considered time scales: On characteristic time scales of the 
deterministic motion only insignificant deviations from the deterministic mo- 
tion typically occur; those trajectories that start close to the boundaries of the 
domains of attraction though are exceptional because they may be influenced 
even by small noise, cross the border of the deterministic domain of attraction 
and, in this way, come close to a "wrong" attractor with finite probability; all 
other trajectories are markedly influenced on much longer time scales only on 
which transitions between the deterministic, locally stable states become likely. 
Hence, these states lose their stability. Nevertheless, for sufficiently weak noise 
the system is found most of the time close to one of the formerly stable states. 
Transitions between these states do occur with certainty even though this hap- 
pens rarely. Therefore such states can be considered as metastable. 

Under the influence of Gaussian white noise the deterministic dynamical 
system (1) becomes a Markov process that is characterized by a Fokker-Planck 
operator of the following form [25, 26] 

^W = -E|:^^(x>*) + E^a.(x,.). (4) 

i i,j ■' 

Wc here will restrict ourselves to periodically driven processes where the drift 
Ki{yi,€) and possibly also the diffusion Dij{x,t) periodically depend on time 
with a common period T. Hence, L{t + T) = L{t). The time evolution of the 
system's probability density function (pdf) p{x,t) is governed by the Fokker- 
Planck equation 

d 

-p{ic,t) = L{t)p{^,t) . (5) 

In the deterministic limit the diffusion matrix vanishes and the drift Ki{x,t) 
approaches the deterministic drift fi{x,t) having the properties discussed in 
Section 1. 

A particular solution of the Fokker-Planck equation is the conditional pdf 
p(x, t\y, s) to find the process at the state x at time t under the condition that 
it was at the state y at time s. It can formally be expressed in terms of the 
Floquet representation in the following way [5, 6, 7, 8, 15] 

p(x, t\y, s) = J2 e''*(*-^Vi(x, t)vi{y, s) , (6) 
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where ijji(x.,t) and ipi{y,s) are Floquet cigcnfunctions and /i^ arc the corre- 
sponding Floquet exponents. They satisfy pairs of mutually adjoint Floquet 
equations reading 

d 

—tpi{x,t) = L{t)ipi{x,t) - iJ.iipi{x,t) , 

g (7) 

-^¥'i(x,t) = L+{t)ipi{x,t) - iJ,iipi{x,t) , 

with natural boundary conditions with respect to the state variable x. Moreover, 
both types of eigenfunctions are periodic in time 

Vi(x,t + T)=Vi(x,t), 
iPi{x,t + T) = ipi{x,t) . 

The Floquet functions tpi{x,t) and (fij{x,t) are mutually orthogonal for eigen- 
values /Uj 7^ fij and can be normalized such that 

dx ipj (x, t)^pi (x, t) = (5i J , (9) 

where 5i,j denotes the Kronecker symbol. The Floquet exponents have real 
parts that are negative or at most zero. 

The representation of the conditional probability in terms of the Floquet 
functions further requires that these functions form a complete set in the sense 
that 

Y,i'i{x,t)^i{Y,t) = 5{x-Y), (10) 

i 

where 5{x) denotes the Dirac 5 function. We note that equations (7), (9) and 
(10) do not uniquely determine the Floquet functions because gauge transfor- 
mations of the form 

i}j{x,t) = gj{t)i)j{x,t) , 

(pj{x,t) = g-\t)ipj{x,t) , ^^^^ 
2m 



with gauge factors 



gj{t) = Cj-e^""^*/^ , Cj € C, Cj ^ (12) 



generate new Floquet eigenfunctions, of. Ref. [27]. Here Z and C denote the 
sets of integer and complex numbers, respectively and i the imaginary unit. 

For the sake of definiteness we assume that the gauge chosen for the Floquet 
representation of the conditional pdf (6) is such that the Floquet exponents 
assume their smallest possible absolute values. The Floquet spectrum consisting 
of these Floquet exponents then contains the value iiq = 0. We assume that 
this Floquet exponent is not degenerate [28] if the diffusion matrix is different 
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from zero. The corresponding eigenfunetion of L~^{i) is constant with respect to 
X and t and can be chosen as ipo{^,t) = 1; the eigenfunetion tpo{yi,t) of L{t) is 
a non-negative and normaUzed function giving the uniquely defined asymptotic 
pdf. Hence, it is the unique solution of the Fokkcr-Planck equation (5) that is 
approached at time t from any initial state in the remote past at s ^ — oo. As 
a Floquet eigenfunetion it is periodic in t. The normalization 

j dxVo(x,i) = l (13) 

follows from eq. (9) together with the fact that (/3o(x, t) = 1. 

For vanishing noise, the diffusion matrix Di j (x, t) vanishes and the backward 
operator becomes a first order partial differential operator Lq (t) = J2i /i(X) t)d/dxi 
with fi{x,t) being the components of the deterministic vector field f(x, gov- 
erning the deterministic motion, cq. (1). For a dynamical system with n coexist- 
ing attractors the characteristic functions of the domains of attraction represent 
n independent periodic solutions of the backward equation —dipo/dt = Lq (t)i^o- 
Each of the solutions is unity on one of the domains of attraction and zero out- 
side. All other periodic solutions are linear combinations of these characteristic 
functions. That means that a deterministic system with n locally stable states 
possesses an n-fold degenerate Floquet eigenvalue /io = 0. As discussed above, 
in the presence of noise, the formerly locally stable states become metastable. 
The n-fold degeneracy of /^o = is lifted, but at sufficiently weak noise there 
remains a group of n Floquet exponents one of which is exactly zero and the oth- 
ers aquire a small negative real part. We call them the slow Floquet exponents. 
For sufficiently small noise this group of slow Floquet exponents stays well sep- 
arated from all other Floquet exponents. For large time lags, the slow Floquet 
exponents and the corresponding Floquet eigenfunctions completely determine 
the conditional pdf which becomes 

ri-l 

p(x,f|y,s) = ^e^'(*-^)Vi(x,t)<^i(y,s) 

i=0 

fort — s ^ T , 

where the sum only runs over the group of n slow Floquet exponents i.e. over 
those exponents with the smallest absolute values. All other Floquet exponents 
are determined by the deterministic time scales all of which are much shorter 
than those given by the slow Floquet exponents. Here r denotes the slowest 
deterministic time scale. 

3.2. Alternative representation of the conditional probability at large times 

In the presence of metastable states the process of moving from a state y at 
time s to a state x at a much later time t may be subdivided into three consecu- 
tive steps that correspond to three contributions to the conditional probability 
p{x,t\y,s): Within the typical relaxation time r, compared to which the con- 
sidered time span t — s is supposed to be very large, the initial state y will 
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be allocated to cither of the mctastablc states (3 with a probability X/3(y, s); 
within the remaining time t — s — TK,t — s the process may visit several other 
metastable states and will be found in the state a at the final time t with a 
probability p{a,t\P, s). Given the final discrete state a, the actual continuous 
states are distributed with a pdf p{^,t\a). For sufficiently small noise the times 
within which the first and the last steps are performed are negligibly short com- 
pared to the total time t — s. Therefore, the initial allocation to a mctastablc 
state a and the final allocation to a continuous state x can be considered as 
instantaneous events. Moreover, all three steps are independent of each other 
and therefore the conditional probability p(x, t\y, s) results as 

p(x, t\y, s) = ^ /9(x, t\a)p{a, t\P, s)x/3(y, s) . (15) 

This particular form of the conditional pdf was derived in the semiadiabatic limit 
[16] which is definded by the regime for which the driving is slow compared to 
the characteristic local relaxation times but not necessarily slow compared to 
the typical transition times between mctastablc states [15]. We claim that this 
particular form of the conditional pdf remains to hold true also beyond the 
semiadiabatic limit, i.e. in situations when the driving period is comparable or 
even faster than the local relaxation times. The rare occurrence of the transi- 
tions between the metastable states is the only condition required for eq. (15) 
to hold. It implies the separation of the times needed to perform the first and 
the third step compared to the much larger time of the second step and justifies 
the independence of these three steps and their respective contributions to the 
conditional probability. Below, we will infer the main properties of these three 
sets of functions p(x, f|a), Xa(x, t) and p{a,t\(3,s) from their according defini- 
tions. 

(i) Each localizing function Xa (x, t) assumes an almost constant value very close 
to unity within the domain of attraction ^^{t) and vanishes outside. Close to 
the border of Vait), the locahzing function Xa{^,i) smoothly interpolates be- 
tween these two values. At each point x all n functions Xa(x, t) exactly add up 
to unity: 

^Xa(x,t) = l. (16) 
a 

(ii) Each a-specific pdf p{x.,t\a) is a strongly peaked function of x about the 
corresponding attractor Aa{t) and rapidly decays away from the attractor. As 
pdf it is normalized to unity 

j dxp{x,t\a) = 1 , (17) 

where the integration extends over the full state space S. Within the respec- 
tive domains of attraction 'Da{t) the a-spccific pdf almost coincides with the 
asymptotic pdf ipoixji) up to a normalizing factor. 

Property (i) of the localizing function allows one to determine the probability 
Pa{t) of finding the metastable state a realized at time t for a given pdf p(x, t) 
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in the following way 

Pa{t) = j (ixxa(x,i)p(x,i) . (18) 

On the other hand, one can assign to a given set of probabiUties Pa{t) a pdf 
Pp(x, t) by decorating the metastable states a with the a-specific pdfs yielding 

Pp(x,t) = 5^p(x,t|a)pc«(t) • (19) 

a 

In order that eqs. (18) and (19) are compatible with each other, i.e. that eq. (18) 
reproduces the prescribed probabilities Pa{t) for p{t) = Pp{t), the localizing 
functions and the a-specific pdfs must form a biorthonormal set of functions, 
i.e. 

dx. x«(x, t)/o(x, t\f3) = Sa.f3 ■ (20) 



For a Fokker-Planck process the time evolution of a pdf p(x, t) is determined by 
the conditional pdf according to 

P{^,t) = j^dy p{x,t\y,s)p{y,s) . (21) 

For large time lags t—s the conditional pdf can be written as in eq. (15). Using 
eqs. (15), (18) and (20) one obtains from eq. (21) for the propagation of the 
probabilities pa {t) 

Pa{t) = Y.p{a,t\0,s)p0{s). (22) 

a,g 

This relation expresses the occupation probabilities of the metastable states at 
a time t in terms of the corresponding probabilities at an earlier time s. Eq. (22) 
hence confirms the interpretation of p{a, s) as the conditional probability of 
the coarse grained process of the metastable, discrete states a = 1 . . .n. 

In order to derive an equation of motion for the probabilities Pa{t) one dif- 
ferentiates both sides of cq. (18) with respect to time, uses the Fokker-Planck 
equation (5), and expresses the pdf by means of eq. (20) in terms of the proba- 
bilities Pfj{t). In this way one obtains 

Mt) = jj^{^-^^^^p{^,t) 

+ Xa(x,t)L(t)p(x,i)} (23) 
= ^K,i}{t)pi}{t) , 

where the time dependent rates kajj{t) arc defined as 



j^dxXa{^,t)L{t)p{x,t\P) . 



(24) 



+ 

/E 
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Eq. (16) implies that the sum over the first index of the rates vanishes, i.e. 
J2a ^oi,i3{t) = 0. Therefore, eq. (23) can be brought into the famihar form of a 
master equation [29] 

= ^",/3(i)P/3(i) - hAt)Pa{t) . (25) 

We expect that for sufficiently low noise the quantities kajj{t) do not become 
negative for a ^ (3 and therefore represent proper rates. A formal proof of 
the positivity though is not available. Negative values of ka,j3{t) though would 
indicate a breakdown of the basic assumption that the long time behavior of 
the process is described by a rate process. 



4. Localizing functions, a-specific pdfs and treinsition rates 

Comparing the two expressions (14) and (15) one finds that the a-spccific 
pdfs p{:x.,t\a) can be expressed as linear combinations of the first n Floquet 
eigenfunctions 'il)i{x.,t) and the localizing functions Xa(x, t) can be written in 
terms of Vi(x, t). This leads to the linear relations 

n-l 

p{x,t\a) = J2Ci,a{t)M^,t), (26) 

i=0 

n-l 

Xa(x,i) = ^Da,i(0^i(x,i), (27) 

i=0 

where Ci^a{t) and Da,i{t) are yet undetermined, time dependent coefficients. 
The orthogonality relations (9), (20) and the linear independence of the first n 
Floquet eigenfunctions imply the following orthogonality relations of the coeffi- 
cients Ci^a{t) and Di^ct{t)'- 

' (28) 
YCiAt)DaAt) = kj ■ 

a 

For i = the normalization of the Floquet function ^o(x, t), see eq. (13), and 
of the a-specific pdfs p{x,t\a), see eq. (17), leads to 

Co,a{t) = 1- (29) 

Next we derive sets of coupled equations of motion for the localizing functions 
and the a-specific pdfs. 
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4.I. Transition rates 

Using the Floquet representation of the a-specific pdfs and localizing func- 
tions, (26) and (27), in combination with the Floquet equations (7) we obtain 
for the rates from eq. (24) 

it)) 

' X ^ (30) 

i 

where we expressed the coefficient matrix Ci^/}{t) as the inverse of D0^i{t) by 
means of eq. (28). Assuming for the moment that the rates ka,/3{t) were known 
we can rewrite eq. (30) as of an equations of motion for the coefficients -Da,i(i) 
and Cj,a(i) reading 

DaAt) = Yl ^c.At)D0At) - DaAt)t^i > (31) 



-CiAt) = Ci,0it)k0At) - f^iCiAt) ■ (32) 



It is interesting to note that these are just the Floquet equations of the 
master equation (25) and, moreover, that the slow Floquet exponents of the 
Fokker-Planck coincide with the Floquet exponents of the master equation. This 
is a consequence of the fact that the master equation specifies the transitions 
between the metastable states, and, therefore, represents the backbone of the 
long time evolution of the Fokker-Planck process. 

With the help of eq. (31) and the Floquet equations (7) the following equa- 
tions of motion for the a-specific pdfs and the localizing functions are obtained 

^p(x, t\a) = L{t)p{^, t\a) - Y ^/3,a(<)p(x, t\f3) , (33) 



-^Xa (X, t) = L+ {t)Xa (X, t)-Y kcc,0 {t)X0 (x, t) . (34) 



These two sets of equations for the functions p(x, t\a) and Xa(x, t) are adjoint 
to each other such that the biorthonormality of the a-specific and the localizing 
functions, sec eq. (20), continues to hold for all times once it holds true at a 
particular instant of time. Eqs. (33) and (34) represent a central result of this 
work. 

The set of coupled equations (33) can be interpreted as the motion of n 
replicas of the original process. Each replica is labeled by one of the attrac- 
tor indices a. The corresponding processes are described by the Fokker-Planck 
equation (5) with additional source and sink terms, X^/j^^^q ^/3 "(^)/'(xj i|a) a-nd 
~'l2i3^a^0A^)p{^'^(^)^ respectively. This means that, say, the a-process dies 
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with probability X!,3/a ^/3.ct(0p(^! ^1/^) '^^'^ instantly resurrects with probabil- 
ity J2(3:^a^0><^i^)Pi^'^\'^) ^^'^^ that the total probability J^, dxp{x, t\a) of each 
replica is conserved for all times. A natural requirement on a process described 
by the set of cqs. (33) is the positivity of the probabilities p(x. t|a). For an ar- 
bitrary choice of the rates ka,/3 {t) this property generally will be violated in the 
course of time. Only for the correct choice of the transition rates the positivity 
is guaranteed to hold. In principle, it is this requirement which determines the 
rates ka,(3{t) on the basis of eq. (33). 

In view of the fact that eqs. (33) and (34) are coupled sets of equations 
not only for the functions p{yi,t\a) and Xc((x, t), respectively, but that in these 
equations also the time dependent rates ka,f3{t) are unknown, it would be very 
difficult to solve these equations exactly. Therefore appropriate approximation 
schemes have to be devised. This will be done in the remaining part of this 
Section. 

4.2. Absorbing boundary approximation: a-specific pdfs 

Assuming the appropriateness of the rate description, i.e. in particular the 
positivity of ka_fj{t) for all a 7^ 0, one can decompose the sum on the right hand 
side of eq. (33) into a sink term ^ X^/j^^^c fc^,a(t) p(x, i|/3) and a source term 
SctT^/j ^fi-a{t)p{'^, t\a). These sink and source terms result from the diagonal and 
non-diagonal parts of the rate matrix {ka,i3{t)), respectively. The sink terms are 
linear combinations of the functions /9(x, f|/3), which are strongly concentrated 
about the positions of the corresponding attractors Ai3{t) with [3 ^ a. 

We approximate these narrow, even though continuously distributed sink 
terms by replacing them with sharp, absorbing states lying on the boundaries 
dBf}{t) of domains Bp{t). Each domain Bfjit) contains the immediate neighbor- 
hood of the attractor Apit) in such a way that the boundary dBf){t) separates 
the corresponding attractor from the remaining state space. Within this ab- 
sorbing boundary approximation we obtain an uncoupled set of equations for 
the a-specific pdfs reading 

—p{y.,t\a) = L{t)p{J^,t\a) + koc{t)p{J^,t\a) , 

forxeS„(t), 
p{x, t\a) = , for all x e dB/s{t) with (3 ^ a , 

where 

fc„(i) = -ka,a{t) = Yl (^^) 

denotes the total decay rate of the state a which is the sum over the individual 

rates from a to all other states /3. The restricted state space '^^{t) is obtained 
from the full state space S by excluding the immediate neighborhoods B/}{t) of 
all metastable states /3 being different from a. Hence, it is defined as 

E„{t) = J:\U0^aBpit). (37) 
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On this restricted state space the function j3{x,t\a) is expected to represent a 
valid approximation of the a-specific pdf p{x,t\a). 

We search for the periodic solution of eq. (35) which can be obtained in the 
following way. First one numerically solves the source free problem 

d 

— p(x, t\a) = L{t)p{x, t\a) , ^^^^ 
p(x, t\a) = 0, for all x G dBfj{t) with l3 ^ a 

with an initial condition that is positive in a small neighborhood of the attractor 
Aa{t) and vanishes everywhere else. Because of the absorbing boundary con- 
ditions at dBi3{t), with /? 7^ a, the auxiliary function p{x,t\a) decays in time, 
i.e. 



Na{t) = / dxp{x,t\a) (39) 
Js„(t) 

is a decreasing function of time. Here the integral is extended over the restricted 
state space 'S.ait) excluding the domains Bi3{t), P ^ a, as defined in eq. (37). 
The normalized function 

p(x, t\a) = p(x, t\a)/Na{t) (40) 
then satisfies the eq. (35) with the total outgoing rate given by 

The such constructed solution p{x,t)/Na{t) approaches a periodic function in 
time on the time scale of the deterministic dynamics, and presents an approxi- 
mation to the a-specific function p{x,t\a). The other rates kp^ait) leaving the 
metastable state a follow from the flux associated with p{x,t\a) through the 
boundaries dB/}{t) 

h,cc{t)= I dS-i{x,t\a), ai^fi, (42) 

JdBp(t) 

where dS denotes the surface element on dB(i{t) pointing towards the metastable 
state Aij{t), and j(x, f|a) the probability current carried by the pdf p{x,t\a). 
Its components read 

3i{x,t\a) = Ki{x,t)p{x,t\a) 

-E^A,Kx,i)p(x,t|a). 

This is a generalization of the well known flux-over-population expression for the 
rate [3, 30, 31, 32]. The stationary flux carrying pdf of the classical flux-over- 
population expression is replaced by the flux carrying time-periodic pdf p{x, t\a) 
which is normalized to one, whence also the population is one. The decisive 
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difference to tlie ciassicai flux-over-population expression lies in the fact that 
in eq. (42) the flux is determined as the probability fiowing per time directly 
into the final metastable state, which because of the surrounding absorbing 
boundary acts as an outlet, rather than through a "saddlcpoint" or "bottleneck" 
on the common part of the separatrices dT>a{t) and d'Dj3{t) of the initial and the 
final metastable state. In the time independent case both expressions coincide 
under the condition that a region containing the final metastable state and 
the bottleneck in question is free of sources [33]. In contrast, in the time- 
periodic case the probability current contains a periodic contribution which 
in general has a nonuniform phase, i.e. the phase depends on the location x. 
Therefore, the instantaneous probability flux through the bottleneck in general 
differs from the flux into the outlet. A large portion of probability flowing 
through the bottleneck, say within the first half of the period may flow back 
during the second half of the period. Only the time averages over one period of 
the probabilities flowing through the bottleneck and into the outlet do coincide. 

4-2.1. a- Floquet functions and rates 

The functions p{x,t\a) which satisfy the Fokker-Planck equation (38) on 
the restricted state space T,a{t) defined in eq. (37) are closely related to the 
Floquet functions ip" (x, t) of the Fokker-Planck operator restricted to {t) 
with absorbing boundaries on the surfaces of the excluded regions Bp{t). These 
a-Floquct functions, as we call them, are the solutions of the corresponding 
Floquet equations which read 

^V?(x,t) = L(t)V?(x,t) - M?^,"(x,i) 

for X e T,a{t) , n = 1,2, 
V-Hx, t) = for X e dB^it), p^a 

Because of the absorbing boundaries at all but one metastable states the Flo- 
quet spectrum consisting of the a-Floquet eigenvalues /if completely lies in the 
complex half plain with negative real part. Wc denote the a-Floquet eigenvalue 
closest to zero by /i" . The absolute value of the real parts of all other a-Floquet 
eigenvalues arc much larger, i.e. <C | for all i ^ 1. In the deterministic 
limit /x" approaches zero, whereas all other a-Floquet eigenvalues stay finite. 
In terms of the a-Floquet eigenfunctions the solution of eq. (38) becomes 

p(x,i|a)=^Cie''?*V'r(x,t), (45) 

i=l 

where Cj are constant coefficients whose values depend on the choice of the initial 

distribution. For times which are large on the deterministic time scale, all terms 
in the sum become negligibly small apart from the first term corresponding to 
Hf. Hence, we obtain 

p(x,t|Q!) oc e'*?* Vf(x,t) , (46) 
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(44) 



and, by proper normalization 



With eq. (41) the total rate ka{t) follows as the negative logarithmic derivative 
of the normalization Jj,^^^ dxtl)"{x,t). It becomes 

fc„(0 = -M? + r„(i), (48) 

where 

r.m = - 

dt 



,{t) = -^Jnl dx^<i{x,t). (49) 

.(*) 



The average of ra{t) over one period vanishes because ra{t) is the derivative of 
a periodic function. Hence, with eq. (48) the a-Floquet eigenvalue /x" is given 

by the negative averaged total rate. 

If one performs the time derivative in eq. (49) one finds 

/^^(,) dx [L(t)V>^(x, t) - (x, t)] 

/s„(t)'^^V'?(x, t) 
IoB,(t) E,, rf^. ai-A.,(x, i)V^?(x, t) (50) 



= E 



= ^fc/3,a(t)+M?. 

In the second equality the time derivative was performed. There, the time 

dependence of the domain T,a{t) does not contribute because the a-Floquct 
function vanishes on the boundary dT,a{t) = y^p^aQB^ii). The time derivative 
of ipf{x,t) was expressed by eq. (44). In the next step the integral involving 
the Fokkcr-Planck operator was written by means of Gauss' theorem in terms 
of surface integrals over the boundary of 'Sa{t). The terms in the sum on /3 are 
the ratios of the probability fluxes through the boundaries dBfs{t) carried by 
the a-Floquet function ijji{x, t), sec eq. (43), and the corresponding populations 
dxipf{x,t). According to the eqs. (42) and (47) the terms in the sum on 
/3 agree with the individual rates kp^a{t)- 



4-3. Absorbing boundary approxim,ations: Localizing functions 

The same type of approximation as for the a-specific pdfs may also be applied 
to the equations of motion for the localizing functions: By neglecting those 
terms on the right hand side of eq. (34) that are proportional to the rates 
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ka.fjit) with (3 ^ a and by introducing absorbing boundary condititions on 
the hypersurfaces dBjs{t)), (3 a we obtain a set of uncoupled equations for 
approximate a-locaUzing functions Xa (x, t) reading 



-^Xa(x,i) 



Xa(x,t) 



WXa(x, t) + ka{t)Xa{^, t) , 

for X e E„ (f) , 
, for ah X e dBnit) with /3 ^ a . 



(51) 



This absorbing boundary approximation is again justified because the rates 
ka,i3(t) are much smaller than the inverse time scales of the deterministic dy- 
namics which govern the motion within the domains of attraction. Moreover it 
is consistent with the above approximation for the a-specific pdfs in the sense 
that the integrals of the products of the a-specific and the respective localizing 
function are independent of time, i.e. 



as follows from eqs. (35) and (51). Note that the time dependence of the in- 
tegration domain T,a{t) docs not contribute because the integrand vanishes at 
the boundary. The biorthonormality of the localizing functions and specific 
pdfs cannot be strictly maintained within this approximation. The deviations 
though are expected to be exponentially small with respect to the noise strength 
because of the small overlap of these functions for different metastable states. 

As in the case of the a-specific functions the total decay rate ka{t) need not 
be known in order to determine the a- localizing functions. Rather one again 
may first determine an auxiliary function Xa(x, t) as the solution of the source 
free equation 



Because of the dissipative nature of the backward operator L~^{t) it is conve- 
nient to integrate this equation backward in time. A forward integration easily 
may run into numerical problems because unavoidable errors would grow ex- 
ponentially in time. As an appropriate final condition for x„(x, to) one may 
choose a function which is constant on the domain of attraction V^ito) and 
zero everywhere else. The solution of this final value problem will approach a 
periodic solution on the time scale of the deterministic dynamics. This asymp- 
totic periodic solution must be normalized at each instant of time by the integral 
of its product with the corresponding a-specific function to yield the required 
approximation of Xa(x, t) 




(52) 



-^Xa(x,t) = L+{t)Xa{x,t) , 

Xa(x, t) = 0, for all x e dB^it) with (3 ^ a . 



(53) 



Xa(x,f) 



Xa(x,t) 
Zait) 



(54) 
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Figure 1: The bistable potential V{x,t), eq. (57), is depicted as a function of the position x 
for different times t = (red, dashed line), t = 0.2T (blue, solid line), and t = OAT (black, 
dotted line) where T denotes the period of the driving and for the driving strength A = 0.1. 

where 

Zait) = dxxa(x,i)p(x,t|a) . (55) 

-'S„(t) 

Using the eqs. (35) and (53) one finds 

*,(0^i|}. (56) 

This relation confirms that the function given by the eqs. (54) and (55) indeed 
is a solution of eq. (51). 

5. Periodically driven Brownian bistable oscillator 

In order to exemplify the theory developed above and to check its consistency 
we consider an overdamped bistable Brownian oscillator driven by an external 
force that varies periodically in time. We choose a bistable quartic potential 
V{x,t) that depends periodically on time, see Fig 1. In conveniently chosen 
dimensionless variables it reads 

V{x,t) =^ -^x'^ + ^x^ - Axsinnt , (57) 

where t is time and x the position of the Brownian particle. The strength of 
the periodic modulation is denoted by A and its frequency by Q. Depending 
on the values of A and f2 the deterministic overdamped dynamics in this time 
dependent potential is either monostable or bistable as displayed in Fig. 2. In 
the present context we are only interested in the bistable region in which the 
deterministic dynamics x = —V'{x,t) possess two stable limit cycles x-i{t) 
and xi{t) and an unstable limit cycle xo{t) forming the separatrix between the 
two attractors, see Fig 3. The diffusion matrix D is taken as constant. The 
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Figure 2: The line dividing the logj^Q f2 — logj^Q A parameter plane into an upper monostable 
and a lower bistable region of the deterministic dynamics x = —V'{x,t) is marked by the 
thick, red solid curve. The blue, thin straight line indicates the value of the forcing strength, 
A'^'^ = 2/{3\/3), below which the potential V{x,t) has two minima for all times t. 




0.25 0.5 0.75 1 

t/T 



Figure 3: The attractors x—i{t), xi(t) and the separatrix xo{t) of the deterministic dynamics 
X = —V'{x,t) for the driving strength A = 0.5 and frequency f2 = 1. 
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Fokker-Planck operator then becomes 



^W = ^^'(^>*) + ^^> (58) 

where V'{x,t) denotes the derivative of the potential with respect to x. The 
corresponding Fokker-Planck and backward equations were numerically solved 
by a collocation method based on a representation of the solution in terms of 
Chebishev polynomials of degree 5 [34]. For all calculations a fixed number 
A'' = 1201 of break-points in the interval [—3, 3] was used. At the ends of 
the interval reflecting boundary conditions were imposed. In the case of the 
forward equation an accuracy of 10~^° led to stable results whereas for the 
backward equations an accuracy of 10~^^ turned out to be necessary in order 
to avoid numerical artefacts. Throughout this paper wc used a fluctuation 
strength given hy D = 1/40. At vanishing driving strength A = the resulting 
bistable symmetric potential then possesses a barrier height per noise energy of 
AV/D = [V{0, 0) - V{1, 0)] /D = 10. 

5.1. Flux- over-population rates 

Wc first numerically determined the time dependent solution p{x, 1 1 — 1) of the 
Fokker-Planck equation (38) on the restricted state space = [— 3,a;i(t)] 

with a reflecting boundary a,t x — —3 and an absorbing boundary at the the 
position of the right attractor Xi {t) and with an initial condition that is sharply 
located at the position of the other attractor (0). After a number n of periods 
T = 27r/0 of the driving frequency had elapsed the remaining population 
N-i{t) was identified as 

/xi{t) 
dxp{x,t\-l), (59) 

see also eq. (39), and the renormalized pdf 

pix,t\-l)^pix,t\-l)/N_iit), (60) 

as well as the rate 

were determined. The number n of transient periods was chosen such that 
ki,-i{t) remained unchanged upon a further increase of n. For different values 
of appropriate numbers n arc collected in Table 1. In Fig. 4 the rates fci__i(i) 
are displayed as functions of time for various driving frequencies. For small 
frequencies the time dependent rate approaches its adiabatic form [16] that is 
given by the inverse mean first time that a process needs to move from x = 
X-i{t) to X = Xi{t) in the frozen potential. The rate then reads [3] 



Jx-i{t) J-3 



(62) 
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Table 1: Number of transient periods 

n n 

1 Too" 

0.5 50 

0.1 10 

0.01 5 

0.001 3 




Figure 4: The rate fci__i(t) following from eq. (61) displays a maximum as a function of 
t/T that becomes lower and shifts towards later times within one period if the frequency fl 
increases. For the frequency O = 10"'^ the rate is indistinguishable from the adiabatic rate 
(62) (black, solid line). The other curves display the rates for Q = 10~^ (blue, dotted line), 0.1 
(red, dash-dotted line), 0.5 (brown, dashed line) and 1 (green, thick dots); in the asymptotic 
limit O — > oo the constant rate A;^^ (thin solid line) given by eq. (63) is approached. In all 
cases the driving strength is A = 0.1 and the noise strength D = 0.025. 



For larger frequencies the maximal value of the rate shrinks and also becomes 
delayed with respect to the driving force. In the limit of high frequencies it 
approches the time independent rate k^^ of a Brownian particle moving in the 
potential V{x,t) = T'^ dt V{x,t) averaged over one period of the driving 
force. For the potential given by eq. (57) the average is symmetric and given 
by V{x,t) — V{x,0). Hence the rate in the limit of high driving frequencies 
coincides with the value of the adiabatic rate at t — 0: 

k^^ = kl'^_^{0) . (63) 

Due to the symmetry of the averaged potential, the rate k^^ also describes the 
opposite transition from the state xi{t) to X-i(t), whence we skipped the index. 

At a fixed frequency the rate fci^_i(i) decreases with decreasing amplitude 
A approaching the time independent value fc'^^, see Fig. 5 

The specific pdf p{x, t\ — 1) given by eq. (60) represents a periodic current 
carrying pdf with an absorbing state at the attractor xi{t). It possesses a single 
maximum the location of which closely follows the deterministic motion of the 
attractor x_i(t), see Fig. 6. The pdf is asymmetric about its maximum with a 
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1 • 10^'^ 



9.75-10- 
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Figure 5; The times at which the rate fci^_i(t) assumes its extrema do hardly depend on the 
ampUtude A. The rate is displayed for various values of A = 0.1 (solid, red), 0.2 (dotted, 
blue), 0.3 (dashdotted, black), and 0.4 (dashed, green); in all cases the frequency is f2 = 10, 
and the noise D = 0.025. Note that for the large amplitude A = 0.4 > A^'^ the deterministic 
attractors x±i{t) are dynamically stabilized, see also Fig. 2. 



breathing width that is wider if the maximum is closer to the position of the 
separatrix xo{t). 

The approximate locahzing function x-ii^Tt) of the left metastable state 
X-i{t) on the restricted state space was obtained from the solution 

X-iix,t) of the backward equation (53) with absorbing boundary condition at 
the right metastable state xi{t). In order to guarantee for sufficient numerical 
stability, the integration of the backward equation has to be performed back- 
ward in time from some to to times t < to- The final function x_i(a;,io) was 
chosen such that it assumes the constant value 1 for all x G [— 3, a:_i(to)] then 
decreases monotonically and reaches zero at the right metastable state. 

After the same number n of transient periods as for the corresponding char- 
acteristic pdf, see Table 1, the normalization integral (55) 

pxi{t) 

Z^i{t)^J dxx~iix,t)p{x,t\-l) (64) 

was determined. The rates fci^_i(t) that follow from the logarithmic derivative 
of Z-i{t), cf. eq. (56), were compared with the rates obtained from eq. (61). 
They are identical within numerical accuracy. 

Finally, the localizing function x-iixTt) was determined by normalizing 
X-i{x,t) with Z-i{t). For an example see Fig. 7. We note that the position 
where the localizing function assumes the value 1/2 coincides with the location 
of the separatrix at the respective time. 

5.2. Floquet approach 

Here we construct the specific pdfs and the localizing functions in terms of 
Floquet eiegenfunctions on the basis of the eqs. (26) and (27). In the present 
case of two metastable states these equations simplify to read 
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-1.5 -1 -0.5 

X 

Figure 6: The specific pdf p{x,t\ — 1) is depicted as a function of the position x for various 
times t = 0.12 T (red, dashed), 0.37 T (blue, soUd), 0.62 T (black, dotted), and 0.87 T (green, 
dashed-dotted) for the driving frequency O = 1, driving amplitude A = 0.1 and noise strength 
D = 0.025. Outside the displayed interval the specific pdf continues to decay. It vanishes at 
the position of the attractor x\{t). The vertical lines indicate the positions of the attractor 
X-\{t) at the respective times. These positions almost coincide with the maxima of the specific 
pdfs at the respective times. 




Figure 7: The localizing function x_i(a::, i) interpolates between the values 1 at the attractor 
X—\{t) and at xi{t). It is displayed at various instants of time, t = 0.12 T (red, dashed), 
0.37T (blue, solid), 0.62r (black, dotted), and 0.87T (green, dash-dotted). The vertical lines 
denote the positions of the separatrix of the deterministic dynamics at the corresponding 
times, see Fig. 3. In the inset a magnification of the center part of the plot marked by a 
rectangle is depicted. It demonstrates that the localizing functions very precisely assume the 
value 1/2 (horizontal line) at the positions of the separatrices indicated by the vertical lines. 
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Table 2: Number of transient periods needed to rea<;h convergence of the Floquet function 
•4>o{x,t) and Floquet exponent ii\ 



9. 


n 


Ml 




1 


10000 


-4.46 10" 


-5 


0.5 


2000 


- 9.46 10" 


-5 


0.1 


1000 


-1.54 10" 


-4 


0.01 


1000 


- 1.58 10- 


-4 


0.001 


100 


-1.58 10- 


-4 



p{x,t\±l) = Mx,t) + C±i{t)i^i{x,t) , (65) 

Here we skipped the first index i of Cj,a(t) since only the values for i = 1 are 
nontrivial in the case of two metastable states. For i = 0, Co,a{t) = 1 always 
holds, see eq. (29). 

To further evaluate these equations (i) the first two Floquet functions of 
the forward and the backward equation and (ii) the coefficients C±i{t) were 
determined numerically The Floquet function ipo{x, t) belonging to the Floquet 
eigenvalue /Uq = is the periodic solution of the Fokker-Planck equation (5), 
(58) with reflecting boundary conditions at x = ±3. As initial condition we 
chose 

^o(.,0)= ^exp(-nx,0)/i.) _ 
J^^dxexp{-Vix,0)/D) 

The Fokker-Planck equation was numerically solved for n periods of the driving 
force. Wc designated this number n in such a way that after subsequent n/10 
periods the Li-norm of the difference of the two solutions was less than 10~^, 
i.e. 

||Vo(a;,l.lnT)-Vo(a;,nT)||i < 10-^ (68) 

where the Li-norm of a function f{x) on the interval [—3,3] is defined by the 
integral of the its absolute value as 

||/(x)||i= f dx\f{x)\. (69) 

The numbers n found in this way are collected in Table 2 for different values 
of the driving frequency. The Floquet function tpi{x,t) and the correspond- 
ing Floquet exponent /ii were obtained from the solution of the Fokker-Planck 
equation (5), (58) with reflecting boundary conditions at a; = ±3 and the initial 
condition 

iJi{x,{))=5{x-X-i{<S)) . (70) 
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After a transient period of duration n T with n given by Table 1 the logarithm 
of the Li-norm of the difference between ipi{x,t) and ipo{x,t) was plotted as a 
function of time for several periods. Its logarithm In 1 1 V'l {x, t) — ipo{x,t)\\\ is the 
superposition of a declining linear and a periodic function of time with period T 
of the driving. The Floquet exponent /ii can be read off from the inclination of 
the linear contribution. The results are presented in Table 2. We note here that 
the method of the a-Floquct fmictions defined on a restricted phase space with 
an absorbing state at, say xi{t), see Section 4.2.1, gave Floquet exponents ii^^ 
which coincide with those based on the full state space up to 4 or 5 digits. The 
same agreement was obtained from the time average of the rates obtained by 
either of the methods described in the previous Section 5.1. Once the Floquet 
exponent fxi is known, the still unnormalized Floquet eigenfunction is obtained 

Mx,t) = e-'''* (M^-,t)-Mx,t) ) . (71) 

The first two Floquet eigenfunctions, which were normalized with respect to the 
Z/i-norm, are displayed in Fig. (8). The Floquet eigenfunction of the backward 
operator belonging to the Floquet exponent /zq = is known to be constant, i.e. 
(po{x, t) = 1. In order to determine the Floquet eigenfunction (pi{x, t) belonging 
to /Ui we solved the backward equation 

-^^^i{x,t) = L+{t)ipi{x,t) (72) 

with the initial condition 

^r(a:,0) = sign(a:).| ^qO • (M - 0.1)^ - 1 \x\toA. ^^^^ 

After a transient time of duration riT with n given in Table 1 all contributions 
from higher Floquet functions have become negligible and if>i{x,t) assumes the 
form 

ipi{x,t) =co + e^'*ciipi{x,t) . (74) 

Knowing the Floquet exponent /Ui we determined the constant such that 
[(fi{x,t) — Co] exp(— /iii) becomes a periodic function of time which is propor- 
tional to the sought-after function ipi{x,t). The normalization of (pi{x,t) is 
chosen such that 

dx (pi{x,t)il)i{x,t) = 1 . (75) 



/: 



The spatial and temporal dependence of (pi{x,t) is depicted in Fig. 9 for the 
same parameter vahics as for the periodic pdf displayed in Fig. 8. 

Once the Floquet functions 'tl)i{x,t) for i = 0, 1 are known the coefficients 
C±i{t) can be determined from the condition that the a-specific pdf p{x,t\a) 
is negligibly small in the vicinity of the other metastable state Xf3{t) {a ^ fi). 
Hence the intergration on both sides of eq. (65) over a small neighborhood 
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Figure 8: The first two Floquet eigenfunctions i/'o(a;,0) (red, solid line) and ipi{x,0) (blue, 
dashed line) of the Fokker-Planck operator (58) of a driven Brownian oscillator in a bistable 
potential (57) for the driving strengths A = 0.1, driving frequency $7 = 1 and noise strength 
D = 2.5 X 10-2 at t = that are displayed in panel (a) are strongly localized in the vicinity 
of the two metastable states at a::±i(0). Both functions are normalized such that their Li- 
norms are one, i.e. | |^/),;(x, t)| |i = dx\ipi(x,t)\ = 1. The two functions almost agree with 
each other up to a change in sign close to the unstable point xo(0). In panel (b), the time 
dependence is indicated for the asymptotic pdf ipo{x,t) for four different times 0.12T (red, 
dashed line), 0.377 (blue, solid fine), 0.62T (black, dotted line) and 0.87T (green, dash-dotted 
line) . 
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Figure 9: The Floquet eigenfunctions (pi(x,t) of the backward operator for the times 0.12T 
(red, dashed Une), 0.37T (blue, soUd Une), 0.62T (black, dotted line) and 0.87T (green, dash- 
dotted line) are almost constant apart from a narrow region about the unstable fixed point 
xo{t). The parameters are with A = 0.1, C = 1 and D = 2.5 X 10~^ the same as in Fig. 8. 



of Xzp{t) gives a negligibly small contribution and thus leads to the following 
expression for the coefficients C±i{t) 

r-x^l(t) + e/2 ^ ihgix, t) 

As an example the coefficient C_i(<) is displayed in Fig. fO for different values of 
the driving frequency. The interval length was chosen as e = O.f . Once the first 
two Floquet eigenfunctions and the coefficients C±i{t) are known, the specific 
pdfs p{x,t\ ± 1) and the localizing functions x±i{^i't) can be calculated and 
compared with the results for p{x,t\ ± 1) and X±i(a;,i), respectively, obtained 
by the flux-over-population method. We here restrict ourselves to a comparison 
for the specific pdf p{x, t\ — f ) for fast driving with $7 = 1. Fig. 11 demonstrates 
the perfect agreement. Only in the immediate vicinity of the metastable state 
a difference becomes visible upon strong magnification. 

Moreover, from the coefficients C±i{t) and the Floquet exponent fii the rate 
and ki^-i(t) can be determined according to eq. (32) which simplifies 
for /ci,_i(i) in the case of two metastable states to 

^^■-^(')= c,{t)^c^,{t) ■ 

A comparison of these rates with those obtained by the reactive flux method is 
presented in Fig. 12 for different values of the driving frequency. A qualitatively 
good agreement is obtained for all frequencies whereby deviations become more 
visible for higher frequencies. 
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Figure 10: The variability of the coefficient C-i(t) within one period T of the driving decreases 
with increasing frequency Q = 10~^ (red, dashed), 10~^ (blue, solid), 10~^ (black, dotted) 
and 1 (green, dash-dotted). The other parameters are with A = 0.1 and D = 2.5 x 10"^ the 
same as in Fig. 8. 




X 



Figure 11: The specific pdf p(x, 0| — 1) was determined by three different methods: As the flux 
carrying periodic pdf p{x,t\ — 1) in the presence of a sharp absorbing boundary at xi{t) (red, 
dashed line), and as a linear combination of the first two Floquet eigenfunctions, see eq. (65), 
with coefficients either determined by eq. (76) (blue, solid line), or from the solutuion of the 
Floquet problem of the master equation (black, dotted line), see the discussion below. Only 
in the magnification displayed in the inset a deviation of the results of these methods becomes 
visible in the vicinity of the metastable state xi(0) ft: 0.98 where p(a;o(0),0| — 1) = 0. We 
expect that these small deviations become even smaller at smaller noise strength. 
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Figure 12: A comparison of the flux-over-population rates (fop rates) (lines) with the Floquet 
rate expressions (F rates) following from eq. (77) (symbols) is presented for frequencies Q = 
0.01 (fop rates: red, solid line; F rates: crosses) and £7 = 0.1 (fop rates: blue, dashed line; F 
rates: circles) in panel (a), and for O = 0.5 (fop rates: red, solid line; F rates: crosses) and 
n = 1 (fop rates: blue, dashed line; F rates: circles) in panel (b). The remaining parameters 
are with A = 0.1, D = 2.5 x 10~^ the same as in Fig. 8. 
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5.3. Decoration 

Finally, we numerically investigated the crucial assumption that after a suf- 
ficiently large transient period the pdf p(x, takes the form of eq. (19), i.e. it is 
determined by the solutions of the master equation (25), Pa{t), which are deco- 
rated by the a-specific pdfs p(x, t\a). As a quantitative measure of the distance 
between the numerically exact solution p{x, t) of the Fokker-Planck equation (5), 
with the Fokker-Planck operator (58), starting at the metastablc state ..c_i(0), 
i.e. with the initial condition (70), and an approximate form pa,{x,t) of the pdf 
we employed the Li-norm (69) of the difference of these functions. The assumed 
asymptotic form 

p^ix,t)=p{x,t\l)pi{t)+pix,t\-l)p-i{t) (78) 

requires the knowledge of the probabilities p±i{t) which was obtained as the 
solution of the master equation 

pi(t) = -fc_i,i(i)pi(<) + fci,_i(t)p_i(t) 
p_i(i) = A;_i,i(t)pi(i) - fci,_i(i)p_i(t) (79) 
Pi(0) = 0, p_i(0) = l, 

where the flux-over-population expressions were taken for the rates, see Section 
5.1. For the specific pdfs we employed three different approximations: First 
we used the current carrying pdfs p{x,t\ ± 1) introduced in Section 5.1. These 
functions were extended onto the full state space [—3, 3] by assigning the value 
zero beyond their respective domains of definition, i.e. we defined 



^ _ / - 3<X< Xl{t) 

PI(X,1\ - I fOTXiit) < x < 3 

0t(x t\l) = ( ° - 3<x<x_i(t) 

Pi[x,T\L) 1^(2.^^11) forx_i(i) <x<3. 



(80) 



As a second and third approximation, in the foUowowing referred to as ap- 
proximation II and III, we used the specific pdfs (65) with the numerically 
determined Floquet functions, see Section 5.2, and determined the coefficients 
C±i{t) in two different ways. The approximation II was obtained by using 
eq. (76) for the coefficients C±i{t). The approximation III is based on the fact 
that these coefficients obey the Floquet equations (32) of the backward master 
equation. We numerically solved these equations under the assumption that the 
rates are given by the flux-ovcr-population expressions. The resulting functions 
c±i{t) then coincide with the sought-after coefficients C±i{t) = qc±i{t) up to 
a common proportionality constant q. Finally this coefficient was determined 
such that the distance between the numerical solution of the Fokker-Planck 
equation and the approximation III, i.e. \\p{x,t) — piu{x,t)\\i, became minimal 
at t = nT with n from Table 1. The coefficients C±i{t) obtained in this way 
are compared with those used in the approximation 11, sec Fig 13. The relative 
deviation between the coefficients C±{t) resulting from the approximations II 



28 



0.996 



1.0002 
1 

0.9998 
0.9996 



1 , 1 

0.25 0.5 


1 

0.75 1 


t/T 




1 1 1 1 1 1 1 


(b) 

■^OOOOOOOOOOOOOOOOOOO 

~x \ 
X X Jk 

X \ / 

X \ M 

X \ / 


ioooooooooS6eoo< 


X \^ ^/x 

X V 
X X 

X X 

^xxxX 
1 . 1 


1 



0.25 



0.5 0.75 

t/T 



Figure 13: The comparison of the approximations II and III for the coefficient C— shows 
perfect agreement for driving frequencies f2 < 0.1, see panel (a) for Q, = 0.1 (method II: 
crosses, method III: solid line). Relatively small but on the scale of the variability apparent 
deviations between the methods become visible for f! = 0.5 (red, method II: crosses, method 
III: solid line) and C = 1 (blue, method II: circles, method III: dashed line) in panel (b). The 
remaining parameters in both panels are with A = 0.1, D = 2.5 X 10~^ the same as in Fig. 8. 
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Figure 14: After a short relaxation time, the decadic logarithm of the L\ distance between 
the numerical solution of the Fokker-Planck equation and the proposed asymptotic form (78) 
reveals a perfect agreement with pni within the expected numerical precision of the solution of 
the Fokker-Planck equation (black, dash-dotted line). In the case of the first method (red, solid 
line) which uses the decoration with the current carrying densities, the absorbing boundary 
conditions at one of the metastable states leads to a larger distance from the asymptotic pdf. 
This also happens with method II (blue, dashed line) which is based on the estimate (76) of 
the coefficients C±\{t) which lacks a rigorous foundation. Yet the observed agreement is very 
good even for rather fast driving with the frequency f2 = 1. The remaining parameters are 
with A = 0.1, D = 2.5 X 10"^ the same as in Fig. 8. 

and III were smaller than 5 x 10^"* in all investigated cases. Clear deviations 
are visible only on the scale of the variability of the coefficients for frequencies 
n > 0.1, see Fig 13. 

The distances between the numerically exact solution of the Fokker-Planck 
equation and the pdfs obtained from the decoration of the metastable states 
according to the three methods described above are displayed in Fig. 14. In 
all cases, after a short initial time, an exponential relaxation sets in until the 
pdfs obtained from method 11 as well as from the decoration with the current 
carrying pdfs saturate at a distance of the order of 2 x 10^''. For method 111 
it does so at the smaller distance of 2 x 10^^. This is a clear indication that 
the asymptotic pdf is indeed of the form of eq. (78). This hence corroborates a 
basic assumption of our work about the structure of the pdf at large times. 

6. Summary 

We investigated the large time stochastic dynamics of periodically driven 
systems with metastable states governed by a Fokker-Planck equation. On time 
scales larger than the typical deterministic time scale this dynamics can be 
completely characterized by the localizing functions, the a-specific pdfs and the 
conditional occupation probabilities of the metastable states. The latter are 
solutions of a Markovian master equation with time-dependent rates. These 
rates can be expressed in terms of the localizing functions and the a-specific 
pdfs, see eq. (24). 
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Using the Floquct representation of the conditional pdf in the large time 
limit we obtained coupled equations of motion for the a-specific densities and 
an adjoint set of equations for the localizing functions. Most interestingly, these 
equations of motion can be interpreted in the spirit of Farkas' [30] and Kramers' 
[31] idea to construct a flux carrying stationary solution by imposing convenient 
sources and sinks. To each a-specific density an a-process can be assigned 
that evolves according to the same dynamical laws as the original process with 
the only difference that it can instantly be translocated in state space. These 
translocations are governed by sinks and sources that cause a sudden death of 
an a-process, say, at a point x and the instant resurrection of the same process 
at a different point y in state space. The sinks are determined by the sum 
of transition rates out of the metastable state a multiplied by those /3 specific 
pdfs corresponding to states that can directly be reached from a. The source is 
given by the total rate to leave state a multiplied by its specific pdf. In this way 
the conservation of probability of each specific pdf is guaranteed. Due to the 
resulting intricate coupling and the dependence on the unknown rates, an exact 
solution is difficult to construct and one must rely on approximate methods to 
solve this set of equations of motion for the a-specific pdfs. 

An efficient way of approximation is based on the fact that at weak noise the 
a-specific pdfs are expected to be strongly localized in the region of the according 
metastable state. This allows one to effectively decouple the equations for the 
a-spccific pdfs (as well as those for the localizing functions) and to calculate a 
current carrying pdf in the presence of sharply absorbing states. The rates of 
all transitions leaving the considered metastable state can then be calculated 
by means of a flux-over-population expression [30, 31, 32]. In contrast to the 
case without time-dependent driving it is important to calculate the probability 
flux flowing directly into the final metastable state. In the time independent 
case this flux is the same through all hypersurfaces in state space separating 
the initial from the final metastable state. In the presence of periodic driving 
the total flux through a hypersurface in general depends both on time and on 
the location of the chosen hypersurface. The proper rate therefore must be 
determined from the probability flux flowing directly into the final metastable 
state. 

We illustrated our theory with the example of a periodically driven bistable 
Brownian oscillator. In contrast to a slowly driven bistable oscillator, at finite 
frequencies bistability extends to larger amplitudes of the driving force. We 
found that the flux-over-population method based on the a-specific pdf with 
an absorbing boundary at the final metastable state requires a much lesser 
computational effort than the direct application of the Floquet approach. In 
the former case the solution of the Fokkcr-Planck equation with the appropriate 
boundary conditions converges on the order of the deterministic time scale, 
whereas for the second method the convergence of the Floquet functions is only 
reached after several transitions between the metastable states have taken place 
on average. 

We note that based on the absorbing boundary approximation the transition 
rates can also be determined by means of numerical simulations of the Langevin 
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equations of the considered Fokkcr-Planck process [17, 18, 19]. 

We finally tested the crucial assumption of our theory saying that the proba- 
bility density resulting as the large time solution of the Fokker-Planck equation 
can be represented as the product of the probabilities of the metastablc states 
decorated by the specific pdfs. The time dependence of the probabilities of the 
metastable states was obtained from the solution of the master equation with 
the numerically determined flux-ovcr-population rates. The specific pdfs ob- 
tained by the absorbing boundary approximation already lead to an excellent 
agreement with the numerically exact solution of the Fokker-Planck equation 
on time scales larger than a few characteristic deterministic times. A more elab- 
orate calculation of the specific pdfs in terms of Floquet eigenfunctions of the 
Fokker-Planck operator led to a further improvement of the agreement by two 
orders of magnitude confirming our assumption. 
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